Discontinuous Shear Thickening in Biological Tissue Rheology

During embryonic morphogenesis, tissues undergo dramatic deformations in order to form functional organs. Similarly, in adult animals, living cells and tissues are continually subjected to forces and deformations. Therefore, the success of embryonic development and the proper maintenance of physiological functions rely on the ability of cells to withstand mechanical stresses as well as their ability to flow in a collective manner. During these events, mechanical perturbations can originate from active processes at the single-cell level, competing with external stresses exerted by surrounding tissues and organs. However, the study of tissue mechanics has been somewhat limited to either the response to external forces or to intrinsic ones. In this work, we use an active vertex model of a 2D confluent tissue to study the interplay of external deformations that are applied globally to a tissue with internal active stresses that arise locally at the cellular level due to cell motility. We elucidate, in particular, the way in which this interplay between globally external and locally internal active driving determines the emergent mechanical properties of the tissue as a whole. For a tissue in the vicinity of a solid-fluid jamming or unjamming transition, we uncover a host of fascinating rheological phenomena, including yielding, shear thinning, continuous shear thickening, and discontinuous shear thickening. These model predictions provide a framework for understanding the recently observed nonlinear rheological behaviors in vivo.


I. INTRODUCTION
During embryonic morphogenesis, biological tissues undergo dramatic deformations in order to form functional organs.The tissues of mature organisms likewise continually suffer stresses and deformations.The success of embryonic development and the maintenance of proper physiological functioning accordingly both depend intimately on a tissue's rheological (deformation and flow) properties [1].On short timescales, tissues can withstand mechanical stresses.Over longer timescales they remodel via cell neighbour exchanges (topological T1 transitions) [2][3][4][5], which thus constitute a key ratelimiting step in important processes such as embryo development, wound healing and cancer metastasis.Recent evidence further suggests that dense confluent tissues, which have no gaps between cells, are poised in the vicinity of a transition between a jammed, solid-like state and an unjammed, fluid-like state [6][7][8][9][10][11][12][13][14].
From a fundamental viewpoint, mechanical stresses can either originate internally within a biological tissue, via spontaneously active processes intrinsic to the cellular level, such as cell contractility [15,16], polarized motility [17] or mitosis [18,19]; or they can be exerted externally, by surrounding tissues and organs [20,21].Recent experiments have shown that cell collectives subjected to externally applied stretching [22][23][24][25] or shear [26] deformations show a strongly non-linear rheological response.Tissues deformed by internally active stresses at the cellular level have likewise been seen to exhibit extreme mechanical phenomena such as fracturing [27].
Perhaps surprisingly, studies of tissue mechanics to date have largely been confined either to the response of tissues to externally imposed stresses or, separately, to phenomena arising from internally active processes.
Crucially, however, most living tissues exist in a state where both forms of driving work together in concert.For example, during Drosophila embryogenesis, polarized actomyosin contractility at the single cell level interacts with external stresses exerted by neighbouring tissues to cause the tissue to flow plastically in convergent extension [21,28].During cancer progression, tumour cell collectives constantly experience mechanical stimuli such as compression and shear stresses from the surrounding extracellular matrix (ECM) [29].At the same time, tumour cells generate actomyosin contractility at the single cell level.This interplay between external microenvironmental stresses and internal motility has been shown to be central to determining whether a cell cluster is jammed or unjammed [12].Recent work on cancer migration also suggests that tumour fluidity depends not only on the single cell invasive potential (akin to our activity) but also on the compressive and shear stresses they experience due to the ECM [30].
With these motivations, in this work, we elucidate in particular the way in which the interplay between globally external and locally internal active driving determines the emergent mechanical properties of the tissue as a whole.Model predictions point towards a framework for understanding the recently observed range of nonlinear rheological behaviors in vivo [27,28,31] and in vitro [23,26].For a tissue in the vicinity of a solidfluid jamming/unjamming transition, we uncover a host of fascinating rheological phenomena, including yielding, shear thinning, continuous shear thickening (CST) and discontinuous shear thickening (DST).
Beyond this context of biological tissues, shear thickening has been the focus of intense recent research in the rheology literature more broadly, due to its widespread occurrence in dense granular materials and suspensions [32][33][34][35].Indeed, simulations [36,37] and experiments [38] on dense suspensions show a large discontinuous increase in viscosity with increasing shear rate, attributed to a cross-over between hydrodynamic and frictional interparticle interactions.For shear rates in this transition region, large stress fluctuations are seen, with an intermittent bimodal switching between low and high viscosity branches of the flow curve [32,37].Associated with this shear thickening transition is the formation of bands of different shear stress stacked with layer normals in the vorticity direction [39].Our prediction of DST in biological tissues suggests that this phenomenon may be present in a broader class of materials than is evident from this existing rheology literature.

II. MODEL
The vertex model that we simulate represents the tightly packed cells of a 2D tissue monolayer as a tiling of n = 1 • • • N polygons, defined by the positions of the polygon vertices [13,[40][41][42][43][44].The vertices of any polygon are joined by edges that form the boundaries with the adjoining cells.Each vertex is shared by three cells and each edge by two cells.
The elastic energy of the vertex model comprises two contributions.The first is set by the deviation of the area A of a cell from a target value A 0 , providing a 2D toy model of 3D cell volume incompressibility.The second contribution is set by the deviation of the cell perimeter P from a target value P 0 .Summing over all cells in the packing, the energy is then The quantity p 0 = P n0 / √ A n0 defines a cell shape factor and is an important control parameter in our study.We set it the same for all cells in any simulation, independent of n.In physical terms, p 0 is commonly attributed to a competition between cell cortical contractility and cellcell adhesion [42,[45][46][47][48], although recent experiments also imply a relationship with cell-substrate traction [49].Cell shape has been shown experimentally to predict jamming behaviour in epithelial tissues [6,50].The elastic constants κ P and κ A set the strength of the perimeter and area interactions, and we choose κ P = 1 as our basic unit of stress.
The elastic forces exerted on the vertices of any cell due to the elastic contributions of that cell are sketched in Fig. 1(a).Consider in this sketch a representative edge of length L, connecting two representative adjacent vertices.The cell that is sketched then contributes to each of these two vertices an equal and opposite tension-like force of magnitude κ P (P − P 0 ), acting tangentially along the edge, inwards along the edge when P > P 0 , and outwards when P < P 0 .The cell shown also contributes to the same two vertices a pressure-like force of magnitude κ A (A − A 0 )L, acting perpendicularly to the edge, in towards the cell when A > A 0 , and outwards when A < A 0 .These expressions are derived in the Appendix.Each vertex in Fig. 1a) additionally belongs to two further cells (not shown) that contribute corresponding elastic forces.The total elastic force ⃗ F j on the jth vertex in the tiling is calculated by summing these contributions from its three shared cells.
In an externally applied simple shear flow of rate γ, with flow direction x and flow-gradient direction y, the position ⃗ r j of the jth vertex in the tiling obeys overdamped dynamics with drag coefficient ζ: with Lees-Edwards periodic boundary conditions [51].
The second term on the right hand side of this equation describes a random motile activity [5,52,53].The magnitude v of this activity is an important control parameter in our study.The direction of the motility of the jth vertex in the tiling is prescribed by the weighted sum of the polarisation vectors ⃗ nij = (cos θ ij , sin θ ij ) of the three cells i j = 1, 2, 3 in contact with that jth vertex.The polarisation angle of each cell in the tiling is initialised at the start of any simulation randomly from in which η n is a random variable with statistics [42] ⟨η The weighting factors w ij in Eqn. 2 ensure that the largest contribution to the polarisation vector of our representative vertex (the jth in the tiling) arises from whichever of its three associated cells i j = 1, 2, 3 has the largest value of the summed lengths of cell edges that contact that vertex.Specifically, we define l ij to be the summed length of the two edges of the i j th cell in contact with vertex j, as shown by the colour coded lines in Fig. 1b), and set consistent with the weighting function used in previous work [5,13].Here L T j is the total length of the three edges in contact with vertex j.Topological T1 cell-cell rearrangement events also intermittently arise, leading to plastic stress relaxation.Specifically, when any cell edge length becomes smaller than a threshold value l T 1 , a T1 event occurs.Prior to a T1, the selected edge is defined by two vertices, one shared between cells αβγ and the other αβδ.The T1 event then replaces these two old vertices with two new ones, shared by cells αγδ and βγδ [13,54].
To initialise an amorphous cellular tiling, we start from a uniform lattice of monodisperse hexagonal cells of cell edge length 1 stacked in √ N rows, each of √ N cells.Target perimeter and area values are then assigned to each cell.To avoid the effects of crystallisation associated with monodisperse packings [55] we use a bidisperse packing in which half the cells have a smaller size and half a larger size.Specifically, we assign these two populations target perimeters in ratio 1 : 1.4 respectively.To maintain a consistent target shape factor (p 0 = P0 √ A0 ) between these two populations, their target areas are set in ratio 1 : 1.4 2 respectively.The overall scale of target area is set such that the target area summed over all cells equals that of the domain size created in the initial uniform hexagonal tiling.The packing is then randomised by implementing cell motility with non-zero v = v prep and D r = D r,prep in the absence of shear for t prep time units then subsequently relaxing the system with zero activity (and zero shear) for t relax time units.Choosing v prep = 4.0, D r,prep = 0.25, t prep = 14.5, t relax = 5.5 was found to produce a random bidisperse initial cellular tiling with fully relaxed cell areas and perimeters.
The equations of motion described above are integrated using the explicit Euler method with time step dt, both during the preparation stage just discussed and the subsequent shearing stage.The timestep dt used was converged to the limit dt → 0, and the system size was converged to the limit N → ∞.
The values, symbols and dimensions for the parameters of the model and shear protocol are listed in Table I.
In what follows, we report the steady state shear stress σ in the tissue as follows.For any individual cell with vertices numbered c = 1 • • • C, a cell level stress is calculated at any time as [56] Here F cα is the α component of the elastic force on vertex c from the elastic contributions of that single cell, and r cβ is the β component of position of vertex c.A is the cell's area.These individual cell stresses, thus defined, are then averaged over all cells in the packing.The packingaveraged shear stress is then averaged over many strain units once a state of statistically steady shear is attained, and furthermore (for each set of model parameters) over three runs with different random number seeds.It is this averaged shear stress that is reported in the results that follow.We have checked it to be robust to changes in system size for N > 100.Fluctuations about the average decrease with increasing N .
In the absence of internal activity (v = 0) and external applied shear ( γ = 0), the vertex model captures a fluid-solid transition at a critical target cell shape p 0 = p * 0 [41,42,48,57], with p * 0 ≈ 3.81 for the bidisperse tiling studied here.For p 0 < p * 0 , cells cannot attain their target shape and the energy barriers to local T1 cell rearrangements are significant: the tissue resists shear, giving a solid phase.For p 0 > p * 0 , cells do achieve their target shape and the energy barriers to rearrangements are small, resulting in a liquid phase that cannot resist shear [48].A nonlinear shear applied quasistatically ( γ → 0) however induces a solidification transition at a critical strain γ c (p 0 ) for p * 0 < p 0 < p * * 0 , with p * * 0 ≈ 4.03 [58].It does so by deforming cells such that they can no longer attain their target shape, eliminating the zero-shear liquid and inducing a solid-like response.The steady state flow curve of shear stress vs shear rate, σ( γ), then displays a yield stress σ Y = lim γ→0 σ( γ) ̸ = 0 for all p 0 < p * * 0 [58].

III. RESULTS
We start by exploring the effects of activity on a sheared tissue, reporting in Fig. 2a) steady state flow curves σ( γ) in the (zero-activity, zero-shear) solid phase, p 0 < p * 0 .At zero activity we see a yield stress, σ Y = lim γ→0 σ( γ) ̸ = 0, indicating a solid-like response with infinite viscosity η = σ/ γ in quasistatic shear γ → 0, consistent with [58].In contrast, at high activity we find liquid-like flow with σ = η γ, in the limit of small shear rate γ → 0. This is termed Newtonian flow behaviour.The viscosity η = η(p 0 , v) is fit by the black dashed lines.
The zero-shear ( γ → 0) viscosity η(p 0 , v) in the (zeroactivity, zero-shear) solid phase, p 0 < p * 0 , thus increases dramatically with decreasing activity v, at fixed target cell shape.In Fig. 3a) we fit this increase to the Vogel-Fulcher-Tamman (VFT) form η ∼ exp [1/(v − v c )] to find the critical activity v = v c (p 0 ) > 0 below which η diverges, at any p 0 < p * 0 .This divergence of the zero shear viscosity at a critical v c (p 0 ) for p 0 < p * 0 may indicate a true yield stress σ Y for all 0 ≤ v < v c , consistent with that at v = 0 [58], although a power law σ ∝ γn with n < 1 is not ruled out for 0 < v < v c .Either way, in the solid phase, p 0 < p * 0 , a critical activity v c (p 0 ) is needed to eliminate solid-like behaviour in favour of Newtonian flow with finite η.This critical v c is plotted vs. p 0 in Fig. 3b), which also shows a colour map of η in the plane of v, p 0 .The zero-shear viscosity is also consistent with the viscosity calculated based on the Green-Kubo relation [59].A linear fit suggests that v c falls to zero at the (zero activity, zero shear) solid-liquid transition p 0 = p * 0 .This intercept is consistent with similar data in linear studies [59].However, the curve shape differs, implying a different mechanism at non-zero activity.Regions of high stress are distributed through the system in d), but formed into system-spanning force chains in e).Target call shape p0 = 3.9, activity v = 0.12.
The flow curves just discussed, for a tissue in its (zeroactivity, zero-shear) solid phase, closely resemble those of complex fluids such as glassy colloidal and jammed athermal soft particle suspensions [60][61][62].These show a yield stress at high packing fraction ϕ and low temperature, analogous to our curves for low activity.Particle suspensions also show low shear Newtonian behaviour at low ϕ and high temperature, analogous to ours at high activity.
We next consider the effect of activity on a sheared tissue in its liquid phase, p 0 > p * 0 See the steady state flow curves in Fig. 2b.In notable contrast to the solid phase, these closely resemble the flow curves of dense frictional suspensions and granular matter [63][64][65].In particular, they show discontinuous shear thickening (DST), in which the shear stress jumps discontinuously with increasing strain rate at high ϕ (in suspensions) or low activity (here).DST then gives way at lower ϕ (in suspensions) or higher activity (here) to continuous shear FIG. 5. Shear rate γthick at the onset of shear thickening as a function of activity for target cell shape p0 = 3.760, 3.800, 3.825, 3.850, 3.875, 3.900, 3.950, 4.00 in curves blue to red leftwards.Coloured straight lines show power law fits to data, γthick ∝ v α .This implies that shear thickening will be present even at very low levels of activity.The dotted black line shows the power α = 2 as a guide to the eye .
thickening (CST), in which the stress still steepens with shear rate, but without jumping.We next explore the origins of DST by examining the spatial-temporal evolution of the stress for different shear rates at fixed activity and p 0 .An example flow curve for v = 0.12, p 0 = 3.90 is shown in Fig. 4a).In Fig. 4b), the stress as a function of strain γ = γt (which is proportional to time t, given constant γ) is plotted for three different shear rates, with line colours corresponding to marker colours in Fig. 4a).The corresponding stress distributions over each of these strain series is shown in Fig. 4b).For low shear rate (blue), the stress fluctuates modestly (in proportional terms) around a low value.Similarly for high shear rate (green) the stress fluctuates modestly around a high value.In contrast, at intermediate shear rate (red), the stress intermittently switches between these low and high stress states to give a bimodal distribution.
This intermittent, bimodal stress evolution at the DST transition is also seen in frictional suspensions, where it is caused by percolating compressive force chains [32-35, 63, 66, 67].Our simulations likewise evidence percolating force chains in this vertex model of biological tissue: Fig. 4c) and d) show representative state snapshots corresponding to the lowest and highest strain rates in Fig. 4a,b), with the thickness of each cell edge proportional to the tensile stress it carries.In the low stress (unthickened) state, the tension is distributed fairly evenly across the system.In contrast, the high stress (thickened) state displays system-spanning force chains.In important contrast to the compressive force chains that form in frictional suspensions, however, we find these stresses to be tensile in nature in tissues.This is consistent with recent computational [10,15,68] and experimental [6,69] studies, which indeed found tensile stresses overwhelmingly to dominate in biological tissues.
To further characterise the regime of shear thickening, we define two shear rates, each via the logarithmic slope of the flow curve, G = d log σ/d log γ. (G = 1 indicates Newtonian behaviour.)First, we define the onset of shear thickening in Fig. 2b) via the shear rate γthick (shown by black squares) at which G first increases above 1 + ϵ, with ϵ = 0.2.Second, we define the reversion to shear thinning at higher strain rates via the shear rate γthin (black triangles) at which G first falls below 1 − δ, with δ = 0.1.At low activity, where DST arises, γthin = γthick , to within the resolution of γ values simulated.Fig. 5) shows γthick as a function of activity v for sev-eral values of the target shape p 0 .For p 0 values comfortably inside the (zero activity, zero shear) fluid phase above p * 0 , we find γthin ≈ γthick ∼ v α at low v.The exponent α decreases with increasing p 0 , with α ≈ 2.0 at p 0 = 4.0.In the fluid phase p 0 > p * 0 , therefore, any level of activity v, however small, is sufficient to restore Newtonian response σ = η γ in quasistatic shear γ → 0, as γ < γthin ≈ γthick ∼ v α .This contrasts notably with the (zero activity, zero shear) solid phase, p 0 < p * 0 , where a finite level of activity v = v c (p 0 ) is needed to give Newtonian behaviour in slow shear, γ → 0.
Having examined the shear rate at which shear thickening (if present) arises in the flow curve for any pairing of p 0 , v values, we consider finally which flow curves indeed show shear thickening.To do so, we define G max to be the maximum of the logarithmic gradient G = d log σ/d log γ across each flow curve and plot in Fig. 6 a colourmap of values of G max that exceed 1 + ϵ with ϵ = 0.2, taking this as the minimal threshold for shear thickening.Values of p 0 , v for which the flow curve does not meet this threshold are shown as white open symbols.As can be seen, very strong shear thickening (large G max ) arises at high p 0 and low v: this is the regime of DST, where the value of G max is limited only by the resolution of γ values simulated.As v increases at fixed p 0 , we see a crossover to more moderate CST before thickening is lost at high v.The black solid line shows a linear fit to the threshold at which thickening is lost.
The apparent loss of shear thickening at fixed v with decreasing p 0 in Fig. 6 is worthy of comment.Towards the left hand edge of the regime of coloured symbols, the shear rate γthick that marks the onset of thickening falls to approach the minimum shear rate that we can feasibly simulate.Were we able to simulate arbitrarily low shear rates, we speculate that the observed regime of thickening would in fact extend leftwards, with γthick → 0 only at the magenta line, consistent with the zero-shear viscosity η being infinite to the left of that line (Fig. 3).We have therefore continued the black solid line leftwards as a dashed line and suggest that the black and magenta lines together delineate the key rheological regimes observed in this work.Representative flow curves for each regime are shown beneath the colourmap in Fig. 6.At higher shear rates, for all p 0 , v, we observe shear thinning arising from T1 cell rearrangement events.This has been seen previously in a vertex model, and derives from an interplay of active fluctuations in vertex length with T1 transitions induced by shear [40].

IV. DISCUSSIONS AND CONCLUSIONS
Our work points towards a framework for understanding the emergent nonlinear mechanics of biological tissue.In particular, we have shown the nonlinear shear rheology of the vertex model to be determined by an intricate interplay between the intrinsic solid-liquid transition that arises at a target cell shape p 0 = p * 0 ≈ 3.81 in the absence of shear or activity [42,70], with the mutually competing effects of a global external shear and local internal cell motility.Indeed, in slow shear, γ → 0, a sufficiently high level of activity always ensures a liquid-like Newtonian response.The path to this liquified state as a function of increasing activity is however markedly different for values of the target cell shape p 0 in the (zero shear, zero activity) solid and liquid phases.In the former, a critical threshold activity level v c (p 0 ) is needed to induce liquefaction.In the latter, any level of activity, however small, ensures Newtonian response in quasistatic shear γ → 0. As the shear rate increases, however, the globally coherent effect of shear exceeds that of locally incoherent activity, inducing a resolidification transition via DST.The shear thickening behavior thus arises from the competition between the accumulation of shear strain due to driving and the dissipation due to cellular activity.On the one hand, the applied shear rate drives the formation of tension networks in the tissue.On the other, the cellular activity acts as a dissipative noisy process that remodels and relaxes the tension network.
We posit that this competition between externally imposed shear and internal cell motility can be characterized via a Peclet number, P e = γτ f [71], in which τ f is the timescale for cell-cell rearrangements due to active motility.In the regime of high Peclet number, P e ≫ 1, γ ≫ τ −1 f , motility is insufficient to affect structural rearrangements caused by the imposed shear.As a result, the mechanical response of the tissue is dominated by the externally applied shear.Because this provides a global driving that acts in a coherent way across the entire tissue, it tends to deform cells away from their target shape, leading to solid-like behaviour.
In contrast, in the regime of low Peclet number, P e ≪ 1, τ −1 f ≫ γ, the mechanics of the tissue is dominated by the active motility.Because this arises internally at the local level of individual cells, lacking any spatial coherence across the tissue, it provides a source of structural rearrangements that tend to counteract the solidifying effect of the globally coherent applied shear just described, resulting in liquid-like response.
We propose that shear thickening occurs at a Peclet number of approximately 1, such that the threshold shear rate for shear thickening is proportional to the inverse of the characteristic time scale τ f .Deriving an expression for τ f is not a simple task.However, we suggest that this timescale for structural rearrangement τ f will be proportional to the tissue's Newtonian viscosity η, defined as the ratio of stress to strain rate in the zero shear rate limit of the flow curve, i.e., at low Peclet number P e ≪ 1.
To explore this idea, we show in Fig 7a) a set of flow curves for a fixed target cell shape p 0 = 3.90, for a range of values of the activity parameter v, now with the shear rate on the horizontal axis rescaled according to γ → γη/v 2 .As can be seen, the location γthick of the shear thickening transition, which is different for different activity values in raw flow curves such as those shown in Fig. 2, now collapses to a single scaled shear rate, 1/ γthick ∝ η/v 2 .This confirms that the inverse shear rate at which thickening occurs in nonlinear rheology, and so our timescale τ f , is indeed proportional to the tissue's zero shear viscosity η.This scaling is further investigated in Fig. 7b), which shows that this relationship 1/ γthick ∝ η/v 2 is approximately obeyed over the full range of values of v and target cell shape p 0 for which both a shear thickening transition and a Newtonian viscosity are indeed seen in the flow curve.
As just described, this rearrangement timescale τ f is important in tissue mechanics because it captures the timescale for structural and stress relaxation driven by internal activity.Via the above scaling argument, we have demonstrated this quantity τ f to be closely related to the tissue viscosity in the limit of zero shear rate.Importantly, this suggests a possible route to accessing the value of τ f experimentally in tissue systems, for example by using magnetically responsive ferrofluid microdroplets to perform quantitative spatiotemporal measurements of mechanical properties in vivo [72,73].
In our current model, the focus has been on the mechanical properties of individual cells and their intercellular interactions, without considering any mechanical role of the cell nucleus.Recent studies have elevated the importance of nuclear compressibility and size as factors that not only govern cell migration and rearrangement but also actively regulate cellular force generation [8,74].In light of these findings, future research should incorporate nuclear mechanics into our existing model, which has already been shown to undergo a density-driven jamming transition [50,75].We anticipate that the model's rheological properties will become increasingly sensitive to nuclear packing density when the size of the nucleus is substantial relative to that of the cell.This could introduce an additional dimension of shear-thickening behavior, similar to the phenomena observed in densely packed particulate systems [63][64][65].
In this work, we have considered a bi-disperse distribution of values of cell size.Looking ahead, it would be valuable to furthermore include phenotypic heterogeneity by incorporating distributions of v, p 0 , κ A and κ P values, grounded in empirical measurements of single-cell properties.Previous research has demonstrated that such heterogeneity can significantly influence tissue rigidity and fluidity [68,76].Consequently, we anticipate that the introduction of mechanical heterogeneity will give rise to intriguing and complex rheological behaviors.
As noted above, DST has been widely observed in dense granular systems with a phenomenology strikingly similar to that reported here for tissues.In each case, a sudden increase in viscosity occurs with increasing shear rate.In each case, this is associated with an intermittent bimodal switching of stress between low and high shear branches, for imposed shear rates in the vicinity of the transition.However, key differences are also notable.In granular systems, DST arises via the development of frictional contacts between particles, leading to the formation of compressive force chains that percolate and bear load across the sample [33,63].In contrast, in tissues we predict DST to arise when the globally coherent effects of an applied shear dominate over the local, spatially incoherent effects of cell motility, leading to the formation of tensile force chains that percolate and bear load.The possibility of vorticity banding associated with DST in tissues should be investigated in future 3D simulations that allow spatial variations in the vorticity direction, to explore the analogy with vorticity banding associated with DST in granular systems [39].
In living tissues, our model predictions can be immediately tested in the convergent extension of the Drosophila germband epithelium [77][78][79].This epithelium experi-ences elongation along the anterior-posterior axis.During this elongation, the germband tissue is subjected to external shearing forces from neighboring structures, such as the ventral furrow, while simultaneously enduring internal forces due to planar polarized contractions driven by myosin II motor activity.This scenario presents a prime example of the interplay between local active forces and global deformations that is central to our theoretical framework.Analyzing the rheological flow curve throughout this process could provide significant insights.With recent technological advances in imaging [80] and force-measurement techniques [81], such analyses are becoming increasingly attainable.
In summary, our study provides a robust framework for understanding the rheological behavior of biological tissues.Considering that nearly all living tissues are subject to a dynamic interplay between local active forces and global deformations, one of the model's most straightforward yet far-reaching predictions is that this interplay can lead to a competition between the timescales of structural relaxation and of external driving forces.This in turn gives rise to diverse rheological responses, including discontinuous shear thickening.Consequently, we argue that these predictions are broadly applicable to a wide array of biological systems.Furthermore, as mounting evidence increasingly suggests that dense tissues operate near a jamming-unjamming transition, our theoretical contributions offer valuable insights into how tissue mechanics is modulated in proximity to these critical states.
In future work, it would be interesting to extend the concepts explored here to understand whether strongly nonlinear mechanical phenomena such as tissue fracture [23,25] and a ductile-to-brittle transition [27] are related to the tissue's ability to shear-thicken.cell's actual area A from its target value A 0 .The second stems from the deviation of the cell's actual perimeter P from its target value P 0 .The resultant force on any given vertex associated with each edge that meets that vertex then likewise comprises separate energy and perimeter contributions, as shown in Fig. 8. (Additional area and perimeter forces also arise associated with the third edge that meets the vertex from neighbouring cells that are not shown.)The area force is shown resolved into components perpendicular to each cell edge.The perimeter force is shown resolved into components parallel to each cell edge.We now derive the magnitude of each of these components.
FIG. 8. Diagram of forces acting on a single vertex from the two edges of one of the cells that meets that vertex.

Area Forces
The cell sketched in Fig. 8 can be separated into an upper triangle and a lower pentagon by the bisecting dotted line shown.With the coordinates of the vertex of interest at (x, y) and the origin (0, 0) defined to coincide with the vertex to its left, the cell area can be represented as the sum of triangle's area and the pentagon's area:

FIG. 1 .
FIG. 1. a) Diagram of elastic forces in the vertex model.Forces tangential to cell edges are proportional to the deviation in cell perimeter and are all of the same magnitude for a single cell.Forces perpendicular to edges are proportional to both the deviation in cell area and the associated edge length.b) Diagram of the jth vertex in the packing (central blue circle), showing the three edges connecting this vertex to its three neighbouring vertices (other blue circles).The cell polarisation vectors ni j of the three adjoining cells ij = 1, 2, 3 are shown as vectors.The associated lengths li j used in the weighted sum to calculate the polarization vector of the vertex are shown by coloured lines.

1 FIG. 2 .
FIG.2.Steady state flow curves of shear stress as a function of shear rate.a) For fixed p0 = 3.65 in the (zero activity, zero shear) solid phase with activity v = 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5 in curves downwards blue to red.For low activity, we find a yield stress in the limit of low strain rate.For high activity, we see Newtonian flow response at low strain rates with shear thinning for higher strain rates.b) For a fixed p0 = 3.90 in the (zero activity, zero shear) liquid phase with activity v = 0.00, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.08, 0.10, 0.12, 0.14, 0.16, 0.20, 0.25, 0.30, 0.35, 0.40 in curves downwards blue to red.With no activity, we find a yield stress in the limit of low strain rate.With modest levels of activity, Newtonian flow response at low strain rate gives way to a discontinuous shear thickening transition with increasing shear rate.Dashed lines fit regimes of constant viscosity, η(p0, v) = σ/ γ.Black squares show γthick and triangles γthin , defined in the main text.

2 FIG. 4 .
FIG. 4. Exploring the discontinuous shear thickening (DST) transition.a) Representative flow curve showing the DST transition.b) Stress as a function of strain γ = γt for the three imposed strain rates denoted by shapes of corresponding colour in a): Blue circles γ = 5.14x10 −5 , red squares γ = 7.36x10 −5 and green triangles γ = 1.1x10 −4 .c) Corresponding probability distributions of the logarithm of the stress.Representative state snapshots at d) γ = 5.14x10 −5 (blue circles) and e) γ = 1.1x10 −4 (green triangles), with the line thickness of any cell edge proportional to the tensile stress it carries.Regions of high stress are distributed through the system in d), but formed into system-spanning force chains in e).Target call shape p0 = 3.9, activity v = 0.12.

FIG. 6 .
FIG. 6. Phase diagram showing different regimes of flow behaviour for different values of v, p0, with a representative flow curve in each regime.Top: Coloured symbols: maximum logarithmic slope Gmax of the flow curve at any v, p0, provided Gmax > 1+ϵ with ϵ = 0.2, designated as the criterion for shear thickening.Open circles have Gmax < 1+ϵ, and no thickening.Black solid line: linear fit to v = v(p0) at which thickening is lost.Dashed line: extrapolation of black solid line left to meet magenta line.Magenta line as in Fig. 3. Panels i-iv) show representative flow curves at the v, p0 values indicated.i) Yield stress, ii) CST, iii) DST, iv) Newtonian.

FIG. 7 .
FIG.7.Data collapse with Peclet number.a) Flow curves for target cell shape p0 = 3.90 and activity v = 0.08, 0.10, 0.12, 0.14, 0.16, 0.20, 0.25, 0.30, 0.35, 0.40 in curves upwards blue to red at left of figure.Compared with flow curves shown in raw form, as for example in Fig.2, the shear rate on the horizontal axis has now been rescaled γ → γ(η/v 2 ) to demonstrate scaling collapse with respect to the location of the shear thickening transition.b) Inverse shear rate at the shear thickening transition 1/ γthick plotted as a function of the scaled viscosity η/v 2 across the full range of values of cell shape p0 and activity v for which a Newtonian regime and a shear thickening transition are observed in the numerically accessible flow curve.The black dashed line shows a linear scaling 1/ γthick ∝ η/v 2 as a guide to the eye.

TABLE I .
Parameters of the model.Dimensions are expressed in terms of modulus (G), time (T) and length (L).